Entrainment and chaos in a pulse-driven Hodgkin-Huxley oscillator 



Kevin K. Lin 

klin@cims . nyu . edu 

January 18, 2006 

Abstract 

The Hodgkin-Huxley model describes action potential generation in certain types of neurons and is a 
standard model for conductance-based, excitable cells. Following the early work of Winfree and Best, 
this paper explores the response of a spontaneously spiking Hodgkin-Huxley neuron model to a periodic 
pulsatile drive. The response as a function of drive period and amplitude is systematically characterized. 
A wide range of qualitatively distinct responses are found, including entrainment to the input pulse train 
and persistent chaos. These observations are consistent with a theory of kicked oscillators developed by 
Qiudong Wang and Lai-Sang Young. In addition to general features predicted by Wang- Young theory, it 
is found that most combinations of drive period and amplitude lead to entrainment instead of chaos. This 
preference for entrainment over chaos is explained by the structure of the Hodgkin-Huxley phase resetting 
curve. 

1 Introduction 

The Hodgkin-Huxley model describes action potential generation in certain types of neurons and is a stan- 
dard model for conductance-based, excitable cells 1 5 18 20 [. There is an extensive literature on the response 
of the Hodgkin-Huxley model to different types of inputs (TlElllIJIilllBlllliZIIlllEllESl, and understand- 
ing how single neurons respond to external forcing continues to be relevant for the study of information 
transmission in neural systems I2T1 1231 . Because neurons typically communicate via pulsatile synaptic 
events, it is natural to investigate the response of the Hodgkin-Huxley model to pulsatile inputs. Early 
studies by Best and Winfree 1 3 38 1 examine the response of a Hodgkin-Huxley model to periodic impulse 
trains, chracterizing in detail the structure of phase singularities and the transition from degree 1 to degree 
phase resetting. However, their work does not systematically address the asymptotic dynamical behavior 
as a function of drive period and amplitude. 1 

This paper studies a spontaneously spiking (i.e. oscillatory) Hodgkin-Huxley neuron model driven 
by periodic, pulsatile input of fixed amplitude and period, and systematically classifies the response as a 
function of drive period and amplitude. It is found that: 

1. In response to periodic pulsatile forcing of fixed amplitude A and period T, a spontaneously spiking 
Hodgkin-Huxley system can exhibit a wide range of distinct behaviors depending on A and T: 

(a) Entrainment: The driven system is stably periodic and its period is a rational multiple of the 
drive period T. 

(b) Transient chaos: The system experiences a transient period of exponential instability before en- 
training to the input. This transient chaos is caused by a Smale horseshoe 1131 . 



Takabe, Aihara, and Matsumoto [32 ] appear to have carried out such a systematic study. But, I was only able to locate an abstract. 
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(c) Chaos: The system becomes fully chaotic: it possesses a positive Lyapunov exponent and a 
mixing attractor (see |39| for a review of these concepts). 

The response of the pulse-driven neuron is approximately To-periodic in the drive period T, where 
Td is the intrinsic period of the unforced Hodgkin-Huxley oscillator. For example, if the pulse-driven 
oscillator is chaotic for some drive amplitude A and drive period T, then it is likely to be chaotic when 
driven by a pulse train of amplitude A with period near T + Tq. 

2. The scenarios enumerated above are prevalent in the sense that they correspond to positive-area sub- 
sets of the drive period-drive amplitude space. Prevalence, together with the approximate periodicity 
stated above, imply that each scenario occurs with positive "probability." (See the discussion of Fig-El 
in [|3|for the precise meaning of probability in this context.) The range of responses and their preva- 
lence are consistent with a theory of nonlinear oscillators developed recently by Qiudong Wang and 
Lai-Sang Young EU EH 

3. While chaotic behavior is readily observable, most combinations of drive period and drive ampli- 
tude lead to entrainment instead of chaos. This preference for entrainment can be explained by the 
structure of the phase resetting curve (see |j4} of the Hodgkin-Huxley system. 

This paper relies heavily on numerical computation and the conceptual framework provided by Win- 
free's theory of biological rhythms |38| and the work of Q. Wang and L.-S. Young on nonlinear oscillators 
1341 1351 . Phase resetting curves, introduced by Winfree, play a particularly important role here. The phase 
resetting curve of a nonlinear oscillator is an interval map which captures the asymptotic response of a non- 
linear oscillator to a single, pulsatile perturbation. Because they are 1-dimensional objects, phase resetting 
curves are often easier to understand than the nonlinear oscillators they represent. They are frequently used 
to infer stable dynamical behavior like phase locking. Wang- Young theory provides a mathematical frame- 
work for using phase resetting curves to infer the existence and prevalence of chaotic behavior. Rather than 
numerically verify the hypotheses of their theorems, we have opted to examine the consequences of the 
theory directly, relying on a combination of numerical simulation and geometric reasoning to characterize 
the specific response of the Hodgkin-Huxley model to a periodic pulsatile drive. 

For the sake of clarity, parameters are selected to ensure that the Hodgkin-Huxley system possesses a 
unique limit cycle and no other attracting invariant set. This corresponds to a repeatedly spiking neuron 
with an unstable rest state. While the scenarios stated above should still hold when the limit cycle coexists 
with other stable invariant sets, this choice simplifies the interpretation of numerical simulations. Other- 
wise, a trajectory may jump out of the basin of the limit cycle, which obscures the mechanism described by 
Wang- Young theory and which Winfree and Best have already investigated thoroughly 131 1381 . 

The rest of this paper is organized as follows: Section |3]briefly reviews the unforced Hodgkin-Huxley 
equations and its properties. Main numerical results are summarized in i|3|and discussed in f|4] Section[5] 
discusses further numerical results, addressing some issues raised in Sections |3] and 0] Section|6| discusses 
possible extensions and generalizations. 

2 Brief review of the Hodgkin-Huxley model 

The Hodgkin-Huxley equations are a system of nonlinear ordinary differential equations 2 which describe 
the way neurons generate spatially and temporally localized electrical pulses \5, 18 20 1. These electrical 

2 This paper does not treat the Hodgkin-Huxley PDEs: spatial dependence is not relevant here. 
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pulses, called action potentials, are the primary way in which neurons transmit information. Action poten- 
tials are triggered by sufficiently large membrane voltages, which can be set up by the influx of ions into 
the cell. A neuron is said to fire or spike when it generates an action potential (Fig.^. The Hodgkin-Huxley 
model describes action potential generation in terms of the membrane voltage and dimensionless gating 
variables which quantify the effective permeability (or conductance) of the membrane for various types of 
ions. 

The original Hodgkin-Huxley equations model action potential generation in the squid giant axon. This 
giant axon contains two types of membrane ion channels. One type of channel is specific to potassium ions, 
the other to sodium ions. The state variables of the model are the membrane voltage V, the activation n of 
the potassium channels, and the activation m and inactivation h of the sodium channels. The equations are 
fT8l 

C^ 1 [-1 - g K n 4 {v - v K ) - g Na m 3 h(v - v Ni ) - g leat (z> - u leak )] 



v 

m = ocJv){\ - m) - p m (v)m 

n = a n (w)(l — n) — (3 n (v)n 

h = a h (v)(l-h)- [3 h (v)h 



(1) 



where 



«ia(*)=v(^)' ^ m (o) = 4exp(i;/18) / 

a» = 0.1W (2+™) , n (v) = 0.125 exp (v/80) , 



(2) 

oc h {v) = 0.07 exp (v/20) , fijv) = l+exp \^y 
^(^) exp(c)— 1' 

Each ion channel consists of independent, identical subunits which must all open to allow ions to pass 
through. The gating variables m, n, and h take value in (0, 1) and represent the fraction of subunits which 
are open. The term n 4 enters into the potassium conductance because potassium channels consist of 4 
identical subunits; analogous structures account for the m 5 h term in the sodium conductance 0. The 
gating variable equations are master equations for continuous-time Markov chains with voltage-dependent 
transition rates a and 13; the Markov chains describe the opening and closing of the corresponding channel 
subunits. The V equation is Kirchoff 's current law. Action potentials are downward voltage spikes and a 
positive I corresponds to an inflow of positively-charged ions. The voltage convention here is that of |18| 
and opposite contemporary usage: the membrane voltage v is defined by 

v = voltage outside — voltage inside. 

Action potentials are generally initiated by perturbations to the membrane voltage. Such perturbations 
may be caused, for instance, by the flow of ions across the cell membrane. Because neurons transmit signals 
through spatially and temporally localized pulses, it is natural to model stimuli as impulses |31|. The 
simplest type of repetitive, pulsatile stimulus to a neuron is a periodic impulse train. This means replacing 
the v equation above by 



-c- 1 



-I - g K n 4 (v - v K ) - g Na m 3 h(v - v Ns ) - g leak (v - v leak ) (3) 



+ A L G(f-fcT), 

where G is a "bump" function such that / G(f) dt = 1. For simplicity, most of this paper uses the choice 
G(t) = 5(t); Section l5~2l discusses the response of the Hodgkin-Huxley system to a pulsatile drive with 

G(t) = I l /t0 ' ° f l ^° • (4) 
0, otherwise 
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Figure 1: The time course for the Hodgkin-Huxley equations at the parameter values 0. The rapid "spike' 
followed by a long "recovery" period is typical of the Hodgkin-Huxley equations. 



Mathematically one can also choose to perturb the gating variables, but such perturbations are not entirely 
natural and are not considered here. 

This paper uses the original Hodgkin-Huxley parameters 1 18 1: 

v Ha = -115 mV, f Na = 120 mQ-Vcm 2 , 

v K = +12 mV, g K = 36 mll^Vcm 2 , ,,-\ 

fie* — —10.613 mV, gteak = 0.3 mCl /cm , 
C = 1 /iF/cm 2 . 

Time is measured in milliseconds and current density in juA/ cm 2 . 

Figure |3 shows a bifurcation diagram for the unforced Hodgkin-Huxley equations. When 1 = 0, the 
neuron maintains a stable rest state, corresponding to the branch of stable fixed points on the left of the 
diagram. A sufficiently large value of I causes a neuron to fire repeatedly, which corresponds to the creation 
of a limit cycle through a saddle-node bifurcation of periodic orbits. Further increasing I destablizes the 
rest state through a subcritical Hopf bifurcation. 

In this paper, the injected current is always set to a value near I « 14, corresponding to a steady ionic 
current which destabilizes the rest state. The phenomena studied here are insensitive to the precise value 
of I as long as it ensures the existence of a stable limit cycle and an unstable fixed point. As explained in 
the Introduction, these properties simplify the interpretation of numerical simulations. For this choice of I, 
the Jacobian of the Hodgkin-Huxley vector field (Eq. E} at the unstable fixed point has two real eigenvalues 
{-4.97815, -0.146991} in the left half plane, and a complex conjugate pair 0.0763367 ± 0.61866/ in the 
right half plane. The fixed point thus has 2-dimensional stable and unstable manifolds. The Lyapunov 
exponents associated with the limit cycle are 0, « —0.20, « —2.0, and « —8.3. Its period is Tq rj 12.944. 
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Figure 2: The bifurcation diagram for the Hodgkin-Huxley equations as the injected current I is varied. The 
line in the middle marks the v coordinate of the rest state. The solid blue part is stable while the dashed red 
part is unstable. Solid black dots near the top and the bottom of the figure are the maximum and minimum 
v values of limit cycles. Empty black circles are the maximum and minimum v values of unstable cycles. 
The fixed point undergoes a subcritical Hopf bifurcation as J increases. This diagram is computed using 
XPPAUT 0. 
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Figure 3: Asymptotic properties of the pulse-driven flow are described by the dynamics of the time-T map 
Ft (see Eq. [5} and its largest Lyapunov exponent A max . Entrainment corresponds to /\ max < 0, and chaos 
corresponds to A max > 0. This figure shows A max as a function of the drive period T, with T ranging from 
Tq 13 (the intrinsic period of the unforced Hodgkin-Huxley system; see jJZJ to 8 • Tq « 101. Left: Kick 
amplitude is A = 5. Right: Kick amplitude is A = 40. Note (i) A max (F + To) « A raax (T); (ii) presence of both 
positive and negative exponents for strong kicks (right), and only zero and negative exponents for weak 
kicks (left); and (iii) the presence of more negative exponents than positive ones. See |]4]for a discussion. 
Lyapunov exponents are estimated by iterating Fj for 1000 steps and tracking the rate of growth of a tangent 
vector. 
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Figure 4: The probability of different response types, as a function of the drive amplitude A. For each 
drive amplitude A, the fraction of drive periods T £ [Tq, 8Tq] for which /\ max (Fj) > 0, etc., is computed by 
sampling from a uniform grid in [Tq, 8Tq\. It is natural to equate these fractions with probabilities because 
the Lyapunov exponents are roughly periodic functions of T (and become more so as X — > oo), as shown 
in FigureOJand explained in Sjl] Empirical definitions: Let A denote the estimated Lyapunov exponent and e 
the estimated standard error. Then "chaos" is defined as A > 3e, "entrainment" A < — 3e, and "rotation" 
\%\ < e/3. 
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3 Main numerical results 



Lyapunov exponents provide a convenient way to characterize the asymptotic dynamics of Eq. [3] Let 
cpt : R 4 — > R 4 denote the flow map generated by the unforced Hodgkin-Huxley equations, T the drive 
period, and A the drive amplitude. The Poincare map 

Fj(v,m,n,h) = (fij(v + A,m,n,h) (6) 

takes a Hodgkin-Huxley state vector (v, m, n,h), applies a pulse of amplitude A to the membrane voltage, 
then evolves it for time T. Iterating the map Fj thus gives a stroboscopic record of the state of our pulse- 
driven Hodgkin-Huxley system before the arrival of each pulse. The long-term dynamical behavior of the 
pulse-driven Hodgkin-Huxley oscillator can be deduced from the asymptotic dynamics of Fj, which is 
characterized by its (largest) Lyapunov exponent /\ max | 13 1: 

A ma< < ■O- Fj has sinks -w- kicked flow is entrained to input 

/\ ma< = ■O- Fj is quasiperiodic •$=>• kicked flow drifts relative to input 
A ma> > ^ Fj chaotic <^> kicked flow is chaotic 

(Sinks refer to stable fixed points and stable periodic orbits.) Note that of the scenarios given in the Intro- 
duction, transient chaos alone does not appear in this list: Lyapunov exponents, being long-time, average 
quantities, cannot detect transient chaos. 

Figure |3 shows the Lyapunov exponents of Fj as a function of T/Tq, where To is the period of the 
Hodgkin-Huxley limit cycle. Different colors correspond to different values of A. The periodicity of the 
response as a function of T is apparent. Because the response type as a function of T is approximately 
identical over each period [nTo, (n + l)To], it makes sense to speak of the probability that a randomly chosen 
drive period T will elicit a particular asymptotic behavior, for example chaos. More precisely, periodicity 
ensures that the fraction p„ of drive periods T in [0, hTq] for which /\ max > converges to a well-defined 
limit as n — > oo. Similar statements hold for /\ max < and /\ max = 0. 

Figure U| shows these probabilities as functions of A. At A = 10, the probability of obtaining a positive 
exponent is roughly 20% and the probability of obtaining a negative exponent is roughly 70%. Thus, if 
one were to pick T randomly out of an interval [NTq, (N + 1)Tq] for large, fixed integer N, the probability 
that /\ max (Fj) > is about 20%. Figure H] shows that when A is small, the most likely type of behavior is 
rotation-like behavior. This possibility becomes less likely as A increases. At the same time, sinks and chaos 
both become more likely, with sinks dominating the scene. One feature of Figure U] specific to the pulse- 
driven Hodgkin-Huxley flow is that when A is large, the system prefers entrainment over chaos in the sense 
that entrainment has higher probability. This preference is more pronounced as A increases. Note that in 
computing Lyapunov exponents numerically, we only have access to finite time information. In principle, 
this means it is virtually impossible to distinguish persistent chaotic behavior from transient chaos caused 
by a "large" horseshoe (but see t i5.lt . 

In all numerical simulations shown in this paper, Eq. is integrated using an adaptive integrator of 
Runge-Kutta-Fehlberg type, with an error tolerance of 10 _6 (in the sup norm) |30|. The largest Lyapunov 
exponent /\ max of Fj is computed in a straightforward manner, by choosing a nonzero unit vector w £ R 4 
and estimating the rate of growth of 1 1 (DFj)"zv\\. The matrix-vector product (DFj) n w is easily computed 
via the variational equations x = H(x),w = DH(x) ■ w for the Hodgkin-Huxley vector field H (DH is the 
Jacobian matrix of H; see 0^1 for details). 
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4 Discussion 



4.1 Response to a single pulse: phase resetting curves 

This section reviews phase resetting curves. See Winfree |38|, Glass and Mackey |11|, and Brown, Moehlis, 
and Holmes |4| for more details and applications, and Guckenheimer and Holmes 1131 for background 
information on dynamical systems theory See |7 8. 9 37 1 for further discussions of phase resetting curves. 

Let (pt '■ R" — > R n be a flow generated by a smooth vector field with a hyperbolic limit cycle y. Such 
a limit cycle represents a stable nonlinear oscillator. The basin of attraction of y is denoted B(y). The 
hyperbolicity of y guarantees that points in B(y) converge to y exponentially fast. (It is convenient to use y 
to refer to both the trajectory y : R — > R" and the invariant point set it defines.) An impulsive perturbation 
("kick") to the nonlinear oscillator can be defined by specifying a kick amplitude A and a kick direction 
K : R n — ► R" and defining a family of kick maps 

K A (x) = x + A -K(x), (7) 

so that kicks send each point x £ R" to K A (x). For what follows, K A should be smooth and satisfy 
K A (B( T )) C B(y). 

The Hodgkin-Huxley system with the value of I given in SjJ] is a nonlinear oscillator whose basin 
B(y) is an open subset of R 4 . The kick map corresponding to an instantaneous voltage spike is simply 
K A (v, m, n, h) = (v + A, m, n, h). As in [J3l it is convenient to introduce the time-T map 

F T = 4> T oK A , (8) 

where o denotes function composition. Iterating Fj gives a stroboscopic record of the system state before the 
arrival of each kick, and thus describes the long-time dynamics of the flow (pt under repeated, F-periodic 
kicks. 

Because the phase dimension n may be large, the dynamics of Fj : R" — > R" may be difficult to analyze. 
Winfree observed that every point near the limit cycle y must converge to y as t — ► oo, so the flow near y is 
dominated by the rotational motion along y. Thus, one can reduce the dimension of the phase space from 
n to 1, at least heuristically. To do this, first define the phase function 9 : y — > [0, To) by fixing a reference 
point xq G y and requiring that (pg( x \(xo) = x for all x £ y. By construction, satisfies j- t (0(y(f))) — 1,0 < 
t < Tq. The function 9 can be extended to a function : B(y) — ► [0, To) by projecting along strong-stable 
manifolds 3 : if y is a point in the basin of y then Q(y) is defined to be 9{x), where x is the unique point such 
that y £ W ss (x). This definition of phase preserves the property that 4r (9((pt{x))) = 1. 

Consider the limit fl2l 

F T = lim F T+nT (10) 

The map Fj is well-defined on the basin of y and retracts the basin onto y, i.e. Fj(x) £ y for all x £ B(y). 
Thus, Fj induces an interval map fj : [0, Fg) — > [0, Fq) which, given the current phase of the system, 
yields the new phase after kicking and evolving the system for time T. That is, fj(9(x)) = 9(Ft(x)) for all 
x £ B(y). 

3 The strong-stable manifold W ss (x) of x £ y is the set 

W ss (x) = {j/eR": lim 0„ To (y)^x). ( 9 ) 

When the vector field generating cpt is smooth, the strong-stable manifolds are (locally) smooth submanifolds of K". The strong-stable 
linear subspace E ss (x) is the tangent space of W ss (x) at x. See IT2IIT3I . 
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Drive amplitude A 


Prob. of sink near plateau 


Prob. of A max < 


5 


41% 


48% 


10 


58% 


62% 


20 


68% 


70% 


30 


76% 


78% 



Table 1: Estimates of the probability of obtaining sinks near the plateau, as a function of A. The data for 
this table is computed by trying about 40 values of T for each choice of A and examining the graph of the 
first return map to the interval [4, 10] (chosen to coincide with the "plateau") and its intersection(s) with 
the diagonal. 

The map fj is the phase resetting curve 4 , or more precisely the finite phase resetting curve (infinitesimal 
phase resetting curves H||7| are not needed here). By construction, it has the property that 

h+s(t) =f T (t)+ 5 (mod T ). (11) 
Thus, the family of maps {fx} is periodic in T. 

Periodicity in drive period T. The approximate periodicity of /\ max (F7-) seen in Figure|3|is easy to under- 
stand heuristically: kicking the oscillator every T seconds and kicking it every T + Tq seconds should yield 
the same asymptotic response because the oscillator simply traverses y at frequency 1/Fo between kicks. 
One can restate this using phase resetting curves: if the drive period T is sufficiently large and 9(xq) = fg, 
then the /j-orbit (fo,/r(fo)//|(fo)/ •■•) should closely follow the phases (8(x ), 8(F T (xo)), 9(Fj(x )), ...) of 
the corresponding Fj-orbit. Since fx+T = /V/ this suggests that A max (Fx+x ) ~ /^(Fr). 

Preference for entrainment. Figure shows phase resetting curves for the Hodgkin-Huxley equations 
for various values of drive period F and drive amplitude A. For sufficientlly small values of A, the phase 
resetting curves are circle diffeomorphisms: either there are sinks (i.e. stable fixed pionts or stable periodic 
orbits), or the map is conjugate to a rotation on a circle and the response of the kicked oscillator drifts 
relative to the periodic drive. As A increases, the graph of fx rather quickly folds over and acquires critical 
points. A striking feature of the graphs in Figure|5]is the "plateau," a phase interval over which fj varies 
very slowly. Another striking feature is the "kink" around 6 ps 9.8. These features are discussed in more 
depth in For now, notice that the plateau provides a simple mechanism for creating sinks: changing the 
kick period T shifts the graph of fj vertically. Whenever the graph intersects the diagonal with a derivative 
\fj\ < 1, then a stable fixed point is created. 

This mechanism can be used to verify the results of Figure|4] compute the graph of the first return map 
of fx to an interval around the plateau, then shift the graph vertically using a number of different values 
of T and estimate the fraction of F's for which fx has a stable fixed point (see Figure |6j. Tabled shows the 
results. For A = 10, the 58% probability of sinks corresponds fairly closely with Figure |4] It is unclear 
whether the ambiguous exponents in Figure |4] really represent positive or negative Lyapunov exponents. 
If a significant fraction of the ambiguous exponents are really negative, then they must come from small 
sinks. 

4 Wang and Young refer to phase resetting curves as singular limits. Phase resetting curves are also sometimes called phase transition 
curves 1111 . 
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Figure 5: The graph of the lift fj of fj, i.e. the unique continuous map K — > K such that fj = fj on [0, To) 
and fj(t + To) = /r(0 (mod To), for the pulse-driven Hodgkin-Huxley equations. Drive amplitudes are 
(a) A = 5, (b) A = 10, and (c) A = 20. The precise value of the drive period T is not so important; varying 
T shifts the graph vertically (Eq. II li . Note that fj has winding number 1 for A = 5 and A = 10, and has 
winding number for A = 20. The numerical data suggests that the degree changes around A ~ 13.589; 
the precise geometric mechanism is not clear. Note that phase ranges from to the intrinsic period To ~ 12.9 
of the Hodgkin-Huxley limit cycle. 
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Figure 6: The first return map Rf T to the interval [4, 10] (chosen to enclose the plateau), for A = 10 and 
T = 17.6. The blue line marks the diagonal. 

Note on numerics. Phase resetting curves are computed here using a variation of the Ermentrout-Kopell 
adjoint method | 4 , SJ. The method is described in the Appendix. A systematic comparison of this method 
to existing methods for computing phase resetting curves is beyond the scope of the present paper and will 
be presented elsewhere. 

4.2 Response to repeated pulses: Wang- Young theory- 
Phase resetting curves provide simple, intuitive explanations for many dynamical properties of pulse- 
driven nonlinear oscillators. For our spiking Hodgkin-Huxley oscillator, explicitly-computed phase reset- 
ting curves show why our pulse-driven neuron prefers entrainment over chaos. In order to infer asymptotic 
behavior, there needs to be a correspondence between the orbits of fj and Fj, and the phase of x E B(y) 
generally does not determine the phase of Fj(x): it may only do so approximately for a finite number of 
iterates. When A max (/Y) < 0, this is enough to show that fj orbits indeed approximate the phases of Fj. 
Inferring chaotic behavior for Fj from fj is far more difficult. Wang- Young theory provides a mathemat- 
ical framework for inferring chaotic behavior using phase resetting curves, and in addition explains why 
chaotic phenomena (and all the other scenarios) is prevalent. 

Shear is an important ingredient of Wang- Young theory. Let y be the limit cycle which represents the 
unforced nonlinear oscillator. Near y, the dynamics follows the periodic, rotational motion on y. Shear 
refers to the presence of an angular velocity gradient around y: the stronger the shear, the sharper the 
angular velocity changes at y. In two dimensions, this means the flow runs much faster on one side of 
y than on the other; in the presence of strong shear, strong stable manifolds tend to become more nearly 
tangent to y. 

Shear and its interaction with kicks is illustrated in a simple model in Figure|7| In the presence of strong 
shear, most ways of kicking the oscillator which take advantage of shear (e.g. kicks which do not take x E y 
too close to the strong-stable manifold W ss (x)) will cause segments of the limit cycle to stretch and fold as 
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(a) (b) 

Figure 7: A simple model (see Eq. 1121 which illustrates the effect of shear. In (a), a large number of initial 
conditions are placed around a limit cycle (the black circle) and a radial kick is applied to each point. The 
blue curve shows the resulting positions of each initial condition. In (b), the kicked points are allowed to 
flow. The red box contains one of the turning points (see ^4.21 . 

they fall back toward y. The phase space stretching caused by shear manifests itself in phase intervals over 
which > 1. This expansion is conducive to chaotic behavior. However, shear also creates region of 
contraction around "turning points," an example of which is highlighted in Figure|7| Such turning points 
correspond to critical points on the phase resetting curve and can easily counteract the expansion needed 
for a positive Lyapunov exponent. The competition between expansion and contraction is the main source 
of difficulty in proving A max (F;r) > 0. 

To infer the existence of parameters for which the time-T map Fj is fully chaotic, Wang and Young use 
results from their previous work on strange attractors with 1 expanding direction and a roughly toroidal 
geometry 1331 1361 . Their main result gives conditions under which there must be T for which A m „(Fj) > 0. 
Furthermore, one can find such "chaotic parameters" near "nice" values of T for which fj has a positive 
Lyapunov expoent. Applying the general theory to kicked oscillators requires checking certain geometric 
conditions. This has been done for a few concrete classes of models 1271 34 1351 . To illustrate the conse- 
quences of the theorems, consider the simple mechanical system 1 34 1 

9(t) + A9(t) = n + A ■ K{6(t)) £ S(t - nT). (12) 

This model was first studied by Zaslavsky, who discovered that this simple system can exhibit fully chaotic 
behavior 1 40 1 . For Eq.^l it can be shown that the full range of scenarios enumerated in the Introduction 
take place and that they are all prevalent. More precisely, Wang and Young prove (see Theorems 1-3 in 1 34 [) 
that: 

1. Invariant curve & weak kicks: When the drive amplitude A is sufficiently small (which is equivalent 
to having a large enough contraction rate X), there exists a simple closed curve y to which all orbits of 
Fj converge and which is invariant under Fj. Moreover, we have the following dichotomy: 

(a) Quasiperiodic attractors: There exist a set of T of positive Lebesgue measure for which Fj is 
topologically conjugate to an irrational rotation. In this case, Fj is uniquely ergodic on f. 
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(b) Gradient-like dynamics: There exists an open set of T such that Fj has a finite number of peri- 
odic sinks and saddles on y, and every orbit converges to one of these periodic orbits. 

2. Gradient-like dynamics without an invariant curve: As A increases (or A decreases), the invariant 
curve y breaks up. Nevertheless, there continues to be an invariant set (no longer a simple closed 
curve) on which gradient-like dynamics persists. 

3. Transient chaos: For even larger A or smaller A, Smale horseshoes (see |13|) will form. Horseshoes 
can coexist with sinks and saddles, creating transient chaos. 

4. Chaos: In the presence of sufficiently strong shear, there exists a positive measure set of drive periods 
T for which Fj is fully chaotic in the sense that it possesses (i) a strange attractor with a positive 
Lyapunov exponent; (ii) at least 1 and at most finitely many ergodic SRB measures 5 with no zero 
Lyapunov exponents; (iii) a central limit theorem; (iv) exponential decay of correlations if a power Fj 
is mixing for some SRB measure V. 

Note that this list of (fairly well-understood) scenarios may not be exhaustive. Other scenarios or combi- 
nations of scenarios are not excluded by the theory. Also, the kicks in Eq.^l ar e purely radial. This is not 
strictly necessary; any kick map which takes advantage of shear will do. See 1 34 1 for precise conditions and 
proofs. 



5 Further results 

5.1 More on the Hodgkin-Huxley phase resetting curve 
Plateau 

The plateau in the phase resetting curve for our pulse-driven Hodgkin-Huxley model (see Figure |HJ cor- 
responds a segment y of the limit cycle y which becomes nearly parallel to a strong-stable manifold after 
receiving a kick. This can be seen by examining the factors which contribute to the derivative fL and which 
can potentially cause fj to become small over a relatively large phase interval. This can be checked by 
writing fj as a composition of other functions and differentiating. 

Let y : R — > R 4 denote the limit cycle trajectory. If we choose y(0) so that 0(y(O)) = 0, then 0(y(f)) = t 
foralH G [0, To), and 

fj = 9 o Fj o y 

= 9 o <p T o K A o y. 

Changing T does not affect fL, so we can set T = 0. Let f — fo- Then f — 9 o K A o y and the chain rule 
gives 

f = (D9oK A oy)-(DK A o r )-f. (13) 
But Ka(v, m, n, h) = (v + A,m,n,h), so its Jacobian DK A is the identity matrix, and for all t S [0, Tg), 

f'(t) = D9(K A (y(t)))-y(t) 

= \D9(K A (y(t)))\ ■ \y(t)\ ■ cos(angle(D0(X A (y(f))),y(f))) (14) 
= \D9(K A (y(t)))\-\y(t)\-sinU(t)), 

5 SRB measures are natural invariant measures for dissipative dynamical systems. They characterize the asymptotic behavior of 
a Lebesgue-positive measure set of initial conditions and have a number of nice mathematical properties. See Young 1391 for an 
introduction. 
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where Z(f) is the angle between y(t) and the strong-stable manifold at K^(y(f)). The last step uses the fact 
that the phase function 9 : B(y) — > [0, Tq) is constant on strong-stable manifolds (see §3). This implies that 
the gradient D8(x) is everywhere orthogonal to W ss (x*), where x* is the unique point in y having the same 
phase as x. The factors in Eq.ll4lthus have simple, geometric meaning: yit) is a point on the limit cycle y 
and \y(t) | is the speed of the limit cycle at that point; is the angle between y(t) and the strong-stable 
manifold at K^(y(f)); and \D9\ measures the rate at which the phase is changing at _K^(y(f)). 

Figure|H]shows /' alongside the 3 factors in Eq.^J The figure shows that the plateau, where fj becomes 
nearly over a long phase interval, coincides with the near-vanishing of The other factors of /' stay 

nearly constant over this interval. Thus, there is a segment y of the limit cycle y, corresponding to the 
phase interval where /.(t) is small, such that K^(y) is nearly tangent to a strong-stable manifold. That 
the segment y, which may be small as a subset of R 4 , corresponds to a large phase interval, is due to the 
relatively slow speed of the limit cycle near y. 

What this argument does not explain is the robustness of this tangency (equivalently, the robustness of 
the plateau) as the drive amplitude A increases (see Fig. This requires a detailed analysis of the geometry 
(in R 4 !) of the strong-stable manifolds (see ij6}. 

Numerical computation of D9. Figure |S] requires the numerical computation of the gradient D9(x) for 
x £ B(y). This can be done as follows: 

Fix x £ B(y) and consider <pj Q (x). Clearly, the limit lim^oo (f>j. (x) = x* exists and has the property 
that x* £ y, 9(x*) = 9(x), and x £ W ss (x*). Set 7r ss (x) = x*. Then n ss projects B(y) onto y and is the identity 
map on y. Furthermore, the nullspace of the Jacobian matrix D7t es (x) of n BS is the tangent to W BB (n ss (x)) at x, 
by construction. 

To compute D9(x), the foregoing discussion suggests that we compute Dip'^ (x) for some large finite n. 
For any finite n, the singular values of Dip j (x) consist of a dominant singular value a\ and 3 nearly zero 
singular values cr 2 , 03, and 04. The a, ; — > as n — > 00 for i = 2,3,4. Denote the left and right singular 
vectors associated with a\ by u and v. It is easy to check that the right eigenvector v is orthogonal to the 
null space of Dn ss (x) and hence tangent to D9. 

This computation requires a relatively accurate estimate of the intrinsic period To of the limit cycle y, 
without which the computation would not converge. This paper adopts the following strategy: instead of 
estimating To just once and reusing its value, solve the system of 24 equations 

x t = H(xi),x 2 = H(x 2 ),J= DH(x 2 ) ■ J (15) 

with initial conditions 

*i(0) G y,x 2 (0) =x r J(0) = Id 4 x4 (16) 

where H is the Hodgkin-Huxley flow field and x is the point at which we would like to evaluate D9. 
Note that Xi,x 2 £ K 4 and / £ R 4x4 . The solution of these equations then give x 2 (t) = <pt{x) and J(t) = 
D4>t(x). The reference trajectory X\ is only used to count the number of periods which have elapsed, and 
the trajectory (x 2 , }) is used to compute Dn ss (x). This procedure works fairly well in practice. 

Kink 

It is natural to ask whether fj (see Fig. [5}, for A — 10 or A = 20, is discontinuous around the kink. 
A discontinuity indicates that there are points in a neighborhood of y which can be kicked outside the 
basin of y. This is not likely the case: Figure shows a magnified view of the phase resetting curve near 
the kink; the graph does not include any numerical interpolation. The figure is obtained by fixing a small 
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Figure 8: The origin of the plateau: in (a), the curves are (i) Black: \f'(t)\; (ii) Blue: |sin(/(i))|; (iii) Red: 
\D9(Kjn(t))\; (iv) Purple: \y(t)\. In (b), the graph of fj near the plateau is shown for reference. Here, 
A = 10. (c) The speed of the Hodgkin-Huxley flow along y. The vertical lines mark the interval [5,9], 
which is part of the plateau. 
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Figure 9: (a) The graph of fj o g^ 1 with drive amplitude A = 10, with the abscissa shown in a new coor- 
dinate system d 1 = g(9) to magnify the region around the "kink." No interpolation is done in this figure: 
only actually computed points are shown, (b) The graph of the coordinate transformation g. The map g is 
generated automatically by the simple adaptive algorithm described in the appendix. 
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Figure 10: A suggestive picture: in (a), a segment of y, starting in the upper left corner of the picture, is 
kicked straight across to the upper right corner. It then follows the flow toward the fixed point for some time 
before spiraling away. The overall direction of motion is top to bottom, (b) Another view of the approach 
to the fixed point. The overall direction of motion here is right to left. The kick amplitude is A = 13.589. 
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Figure 11: Graphs of fj for kick amplitudes A approaching A ait from below (a-c) and above (d-f). The 
horizontal lines mark integral multiples of Tq . 
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parameter 5 > and adaptively refining the grid {9 n } on which the phase resetting curve is evaluated until 
|/r(^n+i) — /t(0»i)| — ^- In Figure^ <5 is set to 0.1. The adaptive procedure (see the Appendix) continues 
to converge for smaller values of 5. 

Figure E| suggests an explanation for the kink: that it is likely caused by a segment of y being kicked 
near the stable manifold of the unstable fixed point. This would cause the segment to wind around the 
stable manifold and eventually spirals away from the fixed point. (Recall that the two unstable eigenvalues 
of the fixed point form a complex conjugate pair.) In the process the kicked segment spreads apart and its 
subsets pick up different amounts of time delays. However, because the Hodgkin-Huxley phase space is 
4-dimensional, FieurellOlcannot give a reliable picture of the dynamics: projecting onto 2 dimensions loses 
too much information. 

The scenario sketched above predicts that there exist a critical kick amplitude A crit at which (y) inter- 
sects the stable manifold of the fixed point. (There may be more than 1 intersection, and more than 1 value 
of A which cause intersections.) As A — > A CIit/ the phase resetting curve should start winding around S 1 
more and more. This can be numerically tested: an estimate of A ait is computed using the Nelder-Meade 
algorithm [30 [ to minimize the closest distance of a trajectory to the fixed point. This yields a critical value 
Aait ~ 13.58953.... When A = A^, exactly, fj should wind around infinitely many times and possess a 
singularity near the location(s) of intersection. For A 7^ A cAt , fj remains smooth, but as A — > A cAt , fj should 
develop a singularity and blow up. See FigurefTTI 

Horseshoes & transient chaos 

Wang- Young theory also guarantees the existence of T's for which Fj exhibits transient chaos, i.e. Fj pos- 
sesses a Smale horseshoe 1131 together with a sink. The coexistence of a horseshoe with a sink has the 
following effect on the dynamics: almost every Fj-orbit would eventually fall into a sink, but an orbit 
which wanders near a horseshoe would dance around unpredictably for a finite number of iterations. Two 
nearby orbits which enter the vicinity of a horseshoe can emerge widely separated, and fall into the sink 
out of phase (unless the sink happens to be a fixed point). In terms of time series data, this kind of behavior 
can be recognized by looking at pairs of trajectories and finding that they chaotically "flutter" about before 
settling down into a steady periodic motion, likely out of phase. 

In contrast to entrainment and chaos, transient chaotic behavior is difficult to observe in the pulse- 
driven Hodgkin-Huxley system. This is because the most likely place to find a horseshoe is near the kink, 
where the expansion so strong that most trajectories escape very quickly. Nevertheless, it is possible to find 
indirect evidence for horseshoes in the pulse-driven Hodgkin-Huxley model. To do so, one looks for an 
interval J C [0, Tq) such that /t(I) gets mapped completely across I at least twice. It is easy, for example, 
to find a "small" horseshoe around the kink in the phase resetting curve. See Figure IT2l The phase interval 
J tells us the rough location of a horseshoe for Fj. 

To go from such an interval I to a horseshoe for the full map Fj, it is necessary to (i) blow up the 
corresponding segment of y to form an open set If C K 4 such that Fj(U) intersects If at least twice, and the 
intersection stretches all the way across If in the unstable direction (along y); and (ii) find invariant cones 
1131 . This can be done in a straightforward manner and is not discussed further here. 

5.2 Miscellaney 

Decay of correlations 

Wang- Young theory predicts that when the dynamics of a pulse-driven oscillator is chaotic and there is a 
unique SRB measure, time correlations (more precisely time autocovariance functions) decay exponentially 
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Figure 12: The phase resetting map and an interval I (marked by the straight lines) which maps across itself 
twice. The abscissa (but not the ordinate) is shown in transformed coordinates. Here, A = 10 and T = 81. 
This is a "small" horseshoe: because the derivative of fj is so large there (on the order of 10 3 ~ 10 4 ), most 
numerically computed orbits escape the horseshoe after a few iterates. 
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Figure 13: Normalized autocorrelation function C m 
indicates that the dynamics is mixing, but the data i 
mixing. 
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(n) for the V variable of the time-T map Fj. This graph 
. insufficient to confirm that the system is exponentially 



fast in time. Figure^lshows the time autocovariance function 

C vv {n) = J o iff) • v A\i - (J v dy\ (17) 



for the voltage variable V (fj. is an ergodic invariant measure). While it clearly decays as n — > oo and thus 
provides evidence that the invariant measure is mixing, the data is not sufficient to confirm that the decay 
is exponential. 



Response to finite-duration pulses 

A natural variation on the numerical experiments of previous sections is to replace instantaneous impulses 
with finite-duration pulses. Heuristically, if the pulse durection fo is less than the fastest of the intrinsic 
timescales of y, the resulting response should be essentially the same as the response to instantaneous 
impulses. With I « 14, these time scales are 12.944 (= the period), 5.1, 0.50, and 0.12 (corresponding to the 
negative Lyapunov exponents). 

FigurelTHsummarizes the numerical results for to = 0.05 (shorter than all time scales), 0.3 (shorter than 
all but one time scale), 2.75 (shorter than all but two fastest time scales), and 9.0 (very slow, not really 
pulsatile in any sense of the word). These graphs should be compared to Figure |I] The pulse amplitude is 
adjusted so that the total amount of charge delivered is the same as an impulse of amplitude A. Interestingly 
enough, the behavior seen earlier are quite robust and disappear only when to = 9. These results suggest 
that the contracting directions do not mix very much over y, and only the slowest contracting time scale 
participates in the production of chaotic behavior. 
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Figure 14: Distribution of Lyapunov exponents for I « 14 with finite-duration pulse of duration (a) fg = 
0.05, (b) to = 0.3, (c) to = 2.75, and (d) to = 9.0. See FigureEJcaption for details. Recall that the Lyapunov 
exponents of the unperturbed limit cycle y are —0.196, —2.01, and —8.31, corresponding to relaxation times 
of 5.1, 0.50, and 0.12. 
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6 Outlook 



The results reported here show that the pulse-driven Hodgkin-Huxley model @ responds to low-frequency 
(relative to the intrinsic period To of the spiking neuron) periodic impulse forcing in a wide range of ways. 
Depending on the drive period and drive amplitude, the response can range from entrainment to fully 
chaotic behavior. This is consistent with the predictions of Wang- Young theory. Furthermore, as shown 
in [01 it is possible to explain some phenomena specific to our pulse-driven Hodgkin-Huxley oscillator in 
terms of special features of the phase resetting curve and provide a partial understanding of the source of 
these features. 

There are some interesting directions for future work: 

Random kick times. The shape of the Hodgkin-Huxley phase resetting curve suggests that if one were 
to drive a Hodgkin-Huxley neuron using a pulse train with random kick times, the resulting random dy- 
namical system can have a negative Lyapunov exponent. This is because the phase resetting curve moves 
up and down from kick to kick, and for any kick time distribution which is sufficiently uniform (e.g. an 
exponential distribution), the probability that the plateau intersects the diagonal is high. The size of the 
plateau suggests that over many iterates, contraction may dominate expansion, leading to a negative Lya- 
punov exponent. A negative Lyapunov exponent implies that two Hodgkin-Huxley neurons, when driven 
by a common pulse train with random kick times, will synchronize. That is, the plateau provides a way to 
create a "random fixed point" [22 \ . These predictions are consistent with preliminary numerical results and 
with a perturbation theory developed by Nakao et. al. for randomly-kicked oscillators in the limit of weak 
kicks 1 26 1. The heuristic, geometric argument sketched above may lead to an extension of their result to the 
regime of strong kicks. 

This synchronization mechanism has also been studied numerically by Doi in the context of a sim- 
ple pieceise linear map |6|. In addition, there is an extensive literature on noise-induced synchrony in 
neural models, including white-noise-driven Hodgkin-Huxley equations 12911411 . Models exhibiting noise- 
induced synchrony provide a concrete framework for exploring neural reliability 12 1 1 1231 . 

Robustness of the phase resetting curve. How robust are the features (plateau, kink) of the Hodgkin- 
Huxley phase resetting curve under perturbations to parameters? How robust is the geometry of the near- 
tangency of kicked segments and strong-stable manifolds responsible for forming the plateau? As the 
range of phenomena predicted by Wang- Young theory may be present in two-dimensional models like the 
Morris-Lecar or FitzHugh-Nagumo, these models may provide a good starting point for exploring these 
questions in lower-dimensional settings. 

A Appendix: Computation of phase resetting curves 

All numerical calculations of phase resetting curves reported in this paper use the algorithm described 
here. It is closely related to an algorithm due to Ermentrout and Kopell |8|. It is presented here for the sake 
of completeness; a systematic comparison with existing methods will appear elsewhere. The algorithm 
computes the phase resetting curve for finite-size perturbations but can be adapted to compute phase resetting 
curves for infinitesimal perturbations. See Appendix A in (H and 1281 1371 1 for more general discussions 
of phase resetting curves, and Ermentrout, Pascal, and Gutkin |9| for a discussion of computing phase 
resetting curves experimentally. 
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The basic idea is to numerically compute the strong-stable linear subspaces along the limit cycle. Then, 
using these linear subspaces as approximations to strong-stable submanifolds, project a kicked point down 
to the limit cycle, where the phase can be estimated. The algorithm really computes strong-stable linear 
subspaces along y and uses these linear subspaces to approximate the phase resetting curve. 

Some preliminaries: the limiting map F T = lim„^oo F T+nTo can be characterized abstractly by the equa- 
tion 

Ft = <Pt ° K A = <p T o 7T ss o K a , (18) 

where 5T ss (x) = y if and only if x E W ss (y), i.e. 7t ee : B(y) — > y maps the basin of y onto the limit cycle 
y along strong-stable manifolds. The abstract notations (and notions) have their uses: the projection 7T SS 
encapsulates the properties of the strong-stable manifolds, e.g. by definition, the strong stable manifold 
W ss (7r ss (x)) passes through x for any x E £>(y), and the nullspace of the Jacobian matrix D7r ss (x) is precisely 
the tangent space of the strong stable manifold W ss (tc ss (x)) at x. 6 

Algorithm (Phase resetting curves via stable subspaces). 

1 . Estimate the period To of the limit cycle y by numerically solving the unforced equations starting with 
a point on or near y. 

2. Discretize y by subdividing the time interval [0, To) into N intervals and computing the correspond- 
ing points Xj E y. Fix an arbitrary reference point xo on y so that each point on y can be assigned a 
unique phase 9 E [0, To). 

3. For each point x, computed in the previous step, compute the Jacobian DH(x,) of the Hodgkin- 
Huxley flow field H at that point. 

4. Using the results of the previous two steps, solve 

x = — H(x), 

£, =j7-(tU>£, (19) 
77 =DH(x) T C 

using the grid points {x, } computed in the previous steps. The x part of the equation above is clearly 
numerically unstable, but that is not a problem because we have already have a numerical represen- 
tation of y. 

The equations above are a variant of the usual method for computing Lyapunov exponents |13|. 
They preserve the length of £,(t), though in practice, it is necessary to rescale £(£) to ensure that 
this constraint is maintained. As t — > oo, £,(t) becomes orthogonal to the strong-stable linear subspace 
E 3S (x(t)) of y. The subspace E ss (x(t)) is tangent to the strong-stable manifold W ss (x(f)) at x(t)7 

6 Note that the commutation relation cpt o 7T SS = 7T SS o cp t expresses the invariance of the strong-stable foliation under cj)f The map 
Ka, in general, has nothing to do with the flow and does not commute with the other maps. 
7 These equations can be generalized to the following: 

x= -H(x) 

h= T}i-{m,li)£.i -Lj<i ((ZUj) + (m,t.j))Zj (20) 
r,, = DH(x) T £; 

If the vectors (£,) form an orthonormal basis at f = 0, then the equations will guarantee that (i;(t)) ar e orthonormal for all t > 0. 
Again, it will be necessary to perform Gram-Schmidt orthogonalizations periodically to maintain this constraint. The vector £j (t), as 
before, converges to a vector orthogonal to E ss (x(t)). So (^2(f)'^3(0/^4(0) span E ss (x(t)). Similarly, the vectors (£3(1), £4(4)) span 
the subspace consisting of the 2 fastest contracting directions, and (£4(4)) spans the fastest contracting direction. 
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5. Using Eq.^Jin combination with the linear subspaces computed in the previous step, we can now 
approximate the phase resetting curve. Start with a point x £ y and compute Of (K^(x)) for increasing 
t. Let to > be the minimum positive time at which Ot (X^(x)) has (i) returned to a small, fixed 
neighborhood of y (in this paper this is chosen to be a neighborhood of distance 10~ 4 around y); and 
(ii) O to (X i 4(x)) lies within one of the pre-computed linear subspaces E ss (x*) for some point x*. Let 0* 
denote the phase of the point x*. Then the new phase of the system is (T + 0* — fp) (mod Tg). 

6. Proceed to the next grid point and repeat. 

7. When the derivative of the phase resetting curve becomes large or infinite, it may be necessary to 
adaptively generate the grid points on which the curve is evaluated. Generally speaking, the grid {x,} 
constructed in Step |3] need not equal the grid {x^} on which the phase resetting curve is evaluated. 
In particular, the grid {x^} can be adaptively chosen to ensure that /«(xj +1 ) — f a {x\) < e, where f a 
denotes the computed phase resetting curve and e is a fixed number, in this paper usually 0.1. This 
adaptive mechanism provides a way to detect discontinuities in fj. 
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